# Dino Christenson & David Glick 
# Replication R script for "Reassessing the Supreme Court: How Decisions & Negativity Bias Affect Legitimacy." Political Research Quarterly.  
# Models 
# June 2018

# wd 
setwd("/prq18_replication/")

# libraries 
library(foreign)

# Load data - PANEL format 
aagmdatawaves <- read.dta("data/CompletePanel_RAVR_Panel_Setup_old.dta")  

#rename and factor these wave variables 
DOMA <- factor(aagmdatawaves$knows_gm_decision_W,
               labels=c("Unaware","Aware"))
VRA <- factor(aagmdatawaves$knows_vra_decision_W, 
              labels=c("Unaware","Aware"))
WAVE <- factor(aagmdatawaves$wave, 
               labels=c("1","2","3","4","5","6"))

#rename issue supports 
gmdoma_sn1 <- aagmdatawaves$gmdoma_index_std_normed_W1
vraaa_sn1 <- aagmdatawaves$vraaa_index_std_normed_W1

#Random Effects via lme4   
library(lme4)
library(pbkrtest)

# for Kenward-Roger p-value approximations 
library(lmerTest)

# individual as random effect

# controls but no interactions (t1.1)
lmer6.noint.did.4W <- lme4::lmer(legit_index_scaled_W ~ 
                                   ideol_distance_W + pid7_W1 +   
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + 
                                   pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 + 
                                   vraaa_sn1 +
                                   gmdoma_sn1 +
                                   WAVE + 
                                   (1|individualid), 
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.noint.did.4W)

# vra policy support interaction (t1.4)
lmer6.vraxw.did.4W <- lme4::lmer(legit_index_scaled_W ~ 
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + 
                                   pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 +
                                   WAVE*vraaa_sn1 + 
                                   (1|individualid), 
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.vraxw.did.4W)

# doma policy support interaction (t1.5)
lmer6.gmxw.did.4W <- lme4::lmer(legit_index_scaled_W ~ 
                                  black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                  sc_info4_index_W1 + media_expo_index_W1 + 
                                  pol_trust_W1 + politicians_disagree_W1 +
                                  gm_family_gay_W1 + 
                                  WAVE*gmdoma_sn1 + 
                                  (1|individualid), 
                                data=aagmdatawaves,
                                REML=T
)
summary(lmer6.gmxw.did.4W)

# party id interactions (t1.2)
lmer6.pidxw.did.4W <- lme4::lmer(legit_index_scaled_W ~ 
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + 
                                   pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 + 
                                   WAVE*pid7_W1 +
                                   (1|individualid), 
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.pidxw.did.4W)

# ideology interactions (t1.3)
lmer6.idlxw.did.4W <- lme4::lmer(legit_index_scaled_W ~ 
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + 
                                   pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 + 
                                   WAVE*ideol_distance_W +
                                   (1|individualid), 
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.idlxw.did.4W) 

# everything interactions (t1.6)
lmer6.allxw.did.4W <- lme4::lmer(legit_index_scaled_W ~ 
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + 
                                   pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 + 
                                   WAVE*pid7_W1 +
                                   WAVE*ideol_distance_W +
                                   WAVE*vraaa_sn1 + 
                                   WAVE*gmdoma_sn1 + 
                                   (1|individualid), 
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.allxw.did.4W) 

# estimate as lme4 (stargazer only reads lme4 not lmerTest)
lmer6.allxw.did.4W.lme4 <- summary(lmer6.allxw.did.4W, ddf = "lme4")
lmer6.allxw.did.4W.lme4

# warning, k-r will take a while after lmerTest... 
# estimates and std errors the same, but t-tests use Satterthwaite approximations to degrees of freedom - a bit more conservative
lmer6.allxw.did.4W.akr <- anova(lmer6.allxw.did.4W, ddf = "Kenward-Roger")
lmer6.allxw.did.4W.kr <- summary(lmer6.allxw.did.4W, ddf = "Kenward-Roger")
lmer6.allxw.did.4W.sw <- summary(lmer6.allxw.did.4W, ddf = "Satterthwaite")
lmer6.allxw.did.4W.kr
lmer6.allxw.did.4W.sw


##### Random Effect Plots #####
par(mfrow=c(1,1))
library(ggplot2)
library(lattice)
# plots qq
qqmath(
  ranef(lmer6.allxw.did.4W, condVar = TRUE), 
  main=F, 
  strip = FALSE)$individualid

# catepillar function approach 
# plot the random effects using caterpillar function
#predictsemer in here
library(AICcmodavg)

Mode <- function(x) {
  ux <- unique(x)
  ux[which.max(tabulate(match(x, ux)))]
}

## re = object of class ranef.mer
ggCaterpillar <- function(re, QQ=TRUE, likeDotplot=TRUE) {
  require(ggplot2)
  f <- function(x) {
    pv   <- attr(x, "postVar")
    cols <- 1:(dim(pv)[1])
    se   <- unlist(lapply(cols, function(i) sqrt(pv[i, i, ])))
    ord  <- unlist(lapply(x, order)) + rep((0:(ncol(x) - 1)) * nrow(x), each=nrow(x))
    pDf  <- data.frame(y=unlist(x)[ord],
                       ci=1.96*se[ord],
                       nQQ=rep(qnorm(ppoints(nrow(x))), ncol(x)),
                       ID=factor(rep(rownames(x), ncol(x))[ord], levels=rownames(x)[ord]),
                       ind=gl(ncol(x), nrow(x), labels=names(x)))
    
    if(QQ) {  ## normal QQ-plot
      p <- ggplot(pDf, aes(nQQ, y))
      p <- p + facet_wrap(~ ind, scales="free")
      p <- p + xlab("Standard normal quantiles") + ylab("Random effect quantiles")
    } else {  ## caterpillar dotplot
      p <- ggplot(pDf, aes(ID, y)) + coord_flip()
      if(likeDotplot) {  ## imitate dotplot() -> same scales for random effects
        p <- p + facet_wrap(~ ind)
      } else {           ## different scales for random effects
        p <- p + facet_grid(ind ~ ., scales="free_y")
      }
      p <- p + xlab("Levels") + ylab("Random effects")
    }
    
    p <- p + geom_hline(yintercept=0)
    p <- p + geom_errorbar(aes(ymin=y-ci, ymax=y+ci), width=0, colour="black")
    # p <- p + geom_point(aes(size=1), colour="blue") 
    p <- p + geom_point(size=0.5, colour="blue") 
    
    # p <- p + theme_bw()	#adding for style but fucks up something
    # p <- p + theme(legend.position="none")
    p <- p + theme(axis.ticks = element_blank(), 
                   axis.text.y = element_blank(), legend.position="none",
                   panel.background = element_rect(fill = "white", colour = NA),
                   panel.border =  element_rect(fill = NA, colour="grey50")
                   # panel.grid.major =  element_line(colour = "grey90", size = 0.2),
                   # panel.grid.minor =  element_line(colour = "grey98", size = 0.5)
                   # panel.margin =   unit(0.25, "lines")
                   # strip.background = element_rect(fill = "grey80", colour = "grey50")
    ) #i'm killing ticks and text to try to make neater for large N small T
    
    return(p)
  }
  lapply(re, f)
}

# ranef qq 
ggcat.qq <- ggCaterpillar(ranef(lmer6.allxw.did.4W, condVar=TRUE))  ## using ggplot2 with qq
ggcat.qq
# just ranef 
ggcat <- ggCaterpillar(ranef(lmer6.allxw.did.4W, condVar=TRUE), QQ=FALSE)  ## this is very nice
ggcat

##### Interaction Plots #####

#plot the interaction with effects package
par(mfrow=c(1,1))
library(effects) 

# (f6.L)
eff.plot.vraxw.did.4W <- plot(
  effect(term="WAVE:vraaa_sn1",
         mod= lmer6.allxw.did.4W, 
         xlevels=list(vraaa_sn1=5), 
         x.var="WAVE",
         se=TRUE, confidence.level=0.95),
  ylim=c(5,10), 
  #xlim=c(2,9),
  ylab="Legitimacy", xlab="Wave", main="",
  key.args=list(x=0.2,y=0.9,corner=c(x=3, y=13)), #positions legend in trellis
  multiline=T)

eff.plot.vraxw.did.4W

# (f6.R)
eff.plot.gmxw.did.4W <- plot(
  effect(term="WAVE:gmdoma_sn1",
         mod= lmer6.allxw.did.4W, 
         xlevels=list(gmdoma_sn1=5), 
         x.var="WAVE",
         se=TRUE, confidence.level=0.95),
  ylim=c(5,10), 
  #xlim=c(2,9),
  ylab="Legitimacy", xlab="Wave", main="",
  key.args=list(x=0.2,y=0.9,corner=c(x=3, y=13)), #positions legend in trellis
  multiline=T)

eff.plot.gmxw.did.4W

# (f5) 
eff.plot.idlxw.did.4W <- plot(
  effect(term="WAVE:ideol_distance_W",
         mod= lmer6.allxw.did.4W, 
         #default.levels=20, 
         xlevels=list(ideol_distance_W=5), 
         x.var="WAVE",
         se=TRUE, confidence.level=0.95),
  ylim=c(5,10), 
  #xlim=c(2,9),
  ylab="Legitimacy", xlab="Wave", main="",
  key.args=list(x=0.2,y=0.9,corner=c(x=3, y=13)), #positions legend in trellis
  multiline=T)

eff.plot.idlxw.did.4W


##### Awareness Models ##### 

# no int (tA2.1)
lmer6.nointawa.did <- lme4::lmer(legit_index_scaled_W ~ 
                                   ideol_distance_W +  pid7_W1 +   
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 + 
                                   vraaa_sn1 + 
                                   gmdoma_sn1 + 
                                   VRA + 
                                   DOMA + 
                                   WAVE + 
                                   (1|individualid),
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.nointawa.did)

# pid int (tA2.2)
lmer6.pidxawa.did <- lme4::lmer(legit_index_scaled_W ~ 
                                  black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                  sc_info4_index_W1 + media_expo_index_W1 + pol_trust_W1 + politicians_disagree_W1 +
                                  gm_family_gay_W1 + 
                                  VRA*pid7_W1 + 
                                  DOMA*pid7_W1 + 
                                  WAVE + 
                                  (1|individualid),
                                data=aagmdatawaves,
                                REML=T
)
summary(lmer6.pidxawa.did)

# idl int (tA2.3)
lmer6.idlxawa.did <- lme4::lmer(legit_index_scaled_W ~ 
                                  black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                  sc_info4_index_W1 + media_expo_index_W1 + pol_trust_W1 + politicians_disagree_W1 +
                                  gm_family_gay_W1 + 
                                  VRA*ideol_distance_W + 
                                  DOMA*ideol_distance_W + 
                                  WAVE + 
                                  (1|individualid),
                                data=aagmdatawaves,
                                REML=T
)
summary(lmer6.idlxawa.did)

# awares * vra supp (tA2.4)
lmer6.vraxawa.did <- lme4::lmer(legit_index_scaled_W ~ 
                                  black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                  sc_info4_index_W1 + media_expo_index_W1 + pol_trust_W1 + politicians_disagree_W1 +
                                  gm_family_gay_W1 +
                                  VRA*vraaa_sn1 + 
                                  DOMA*vraaa_sn1 + 
                                  WAVE + 
                                  (1|individualid),
                                data=aagmdatawaves,
                                REML=T
)
summary(lmer6.vraxawa.did)

# awares * doma supp (tA2.5)
lmer6.gmxawa.did <- lme4::lmer(legit_index_scaled_W ~ 
                                 black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                 sc_info4_index_W1 + media_expo_index_W1 + pol_trust_W1 + politicians_disagree_W1 +
                                 gm_family_gay_W1 +
                                 VRA*gmdoma_sn1 + 
                                 DOMA*gmdoma_sn1 + 
                                 WAVE + 
                                 (1|individualid),
                               data=aagmdatawaves,
                               REML=T
)
summary(lmer6.gmxawa.did)

# both vra*aa and doma*gm ==all waves with full controls too 
lmer6.suppxawa.did <- lme4::lmer(legit_index_scaled_W ~ 
                                   ideol_distance_W + pid7_W1 +   
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 + 
                                   VRA*vraaa_sn1 + 
                                   DOMA*gmdoma_sn1 + 
                                   WAVE + 
                                   (1|individualid),
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.suppxawa.did)

# all int (tA2.6)
lmer6.allxawa.did <- lme4::lmer(legit_index_scaled_W ~ 
                                  black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                  sc_info4_index_W1 + media_expo_index_W1 + pol_trust_W1 + politicians_disagree_W1 +
                                  gm_family_gay_W1 +
                                  VRA*pid7_W1 + 
                                  DOMA*pid7_W1 + 
                                  VRA*ideol_distance_W + 
                                  DOMA*ideol_distance_W + 
                                  VRA*vraaa_sn1 + 
                                  DOMA*vraaa_sn1 +
                                  VRA*gmdoma_sn1 + 
                                  DOMA*gmdoma_sn1 +
                                  WAVE + 
                                  (1|individualid),
                                data=aagmdatawaves,
                                REML=T
)
summary(lmer6.allxawa.did)


# ranef qq 
ggcat.awa.qq <- ggCaterpillar(ranef(lmer6.allxawa.did, condVar=TRUE))  ## using ggplot2 with qq
ggcat.awa.qq
# just ranef 
ggcat.awa <- ggCaterpillar(ranef(lmer6.allxawa.did, condVar=TRUE), QQ=FALSE)  ## this is very nice
ggcat.awa


##### Interaction Plots for Awareness #####

#plot the interaction with effects package
par(mfrow=c(1,1))
# library(effects) 

# (fA4.C)
eff.plot.vraxawa.did <- plot(
  effect(term="VRA:vraaa_sn1",
         mod= lmer6.allxawa.did, 
         #default.levels=20, 
         xlevels=list(VRA=2), 
         x.var="vraaa_sn1",
         se=TRUE, confidence.level=0.95),
  ylim=c(6,10), 
  #xlim=c(2,9),
  ylab="Legitimacy", xlab="VRA Support", main="",
  key.args=list(x=0.2,y=0.9,corner=c(x=3, y=13)), #positions legend in trellis
  multiline=T)

eff.plot.vraxawa.did

# (fA4.R)
eff.plot.gmxawa.did <- plot(
  effect(term="DOMA:gmdoma_sn1",
         mod= lmer6.allxawa.did, 
         #default.levels=20, 
         xlevels=list(DOMA=2), 
         x.var="gmdoma_sn1",
         se=TRUE, confidence.level=0.95),
  ylim=c(5,10), 
  #xlim=c(2,9),
  ylab="Legitimacy", xlab="GM Support", main="",
  key.args=list(x=0.2,y=0.9,corner=c(x=3, y=13)), #positions legend in trellis
  multiline=T)

eff.plot.gmxawa.did 

# (fA4.L)
eff.plot.idlxawa.did <- plot(
  effect(term="VRA:ideol_distance_W",
         mod= lmer6.allxawa.did, 
         #default.levels=20, 
         xlevels=list(VRA=2), 
         x.var="ideol_distance_W",
         se=TRUE, confidence.level=0.95),
  ylim=c(5,10), 
  #xlim=c(2,9),
  ylab="Legitimacy", xlab="VRA Support", main="",
  key.args=list(x=0.2,y=0.9,corner=c(x=3, y=13)), #positions legend in trellis
  multiline=T)

eff.plot.idlxawa.did

##### Single Variable Issue Support Model #####

# everything interactions (tA3.6)
lmer6.allxw.did.iss <- lme4::lmer(legit_index_scaled_W ~ 
                                    black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                    sc_info4_index_W1 + media_expo_index_W1 + 
                                    pol_trust_W1 + politicians_disagree_W1 +
                                    gm_family_gay_W1 + 
                                    WAVE*pid7_W1 +
                                    WAVE*ideol_distance_W +
                                    WAVE*vra_states_W1 + 
                                    WAVE*gm_hetero_likert_W1 + 
                                    (1|individualid), 
                                  data=aagmdatawaves,
                                  REML=T
)
summary(lmer6.allxw.did.iss) 

# controls but no interactions (tA3.1)
lmer6.noint.did.iss <- lme4::lmer(legit_index_scaled_W ~ 
                                    ideol_distance_W + pid7_W1 +   
                                    black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                    sc_info4_index_W1 + media_expo_index_W1 + 
                                    pol_trust_W1 + politicians_disagree_W1 +
                                    gm_family_gay_W1 + 
                                    gm_hetero_likert_W1 + 
                                    vra_states_W1 +
                                    WAVE + 
                                    (1|individualid), 
                                  data=aagmdatawaves,
                                  REML=T
)
summary(lmer6.noint.did.iss)

# vra policy support interactions (tA3.4) 
lmer6.vraxw.did.iss <- lme4::lmer(legit_index_scaled_W ~ 
                                    black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                    sc_info4_index_W1 + media_expo_index_W1 + 
                                    pol_trust_W1 + politicians_disagree_W1 +
                                    gm_family_gay_W1 +
                                    WAVE*vra_states_W1 + 
                                    (1|individualid), 
                                  data=aagmdatawaves,
                                  REML=T
)
summary(lmer6.vraxw.did.iss)

# doma policy support interactions (tA3.5)
lmer6.gmxw.did.iss <- lme4::lmer(legit_index_scaled_W ~ 
                                   black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                   sc_info4_index_W1 + media_expo_index_W1 + 
                                   pol_trust_W1 + politicians_disagree_W1 +
                                   gm_family_gay_W1 + 
                                   WAVE*gm_hetero_likert_W1 + 
                                   (1|individualid), 
                                 data=aagmdatawaves,
                                 REML=T
)
summary(lmer6.gmxw.did.iss)

# party id interactions (tA3.2) 
lmer6.pidxw.did.iss <- lme4::lmer(legit_index_scaled_W ~ 
                                    black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                    sc_info4_index_W1 + media_expo_index_W1 + 
                                    pol_trust_W1 + politicians_disagree_W1 +
                                    gm_family_gay_W1 + 
                                    WAVE*pid7_W1 +
                                    (1|individualid), 
                                  data=aagmdatawaves,
                                  REML=T
)
summary(lmer6.pidxw.did.iss)

# ideology interactions (tA3.3)
lmer6.idlxw.did.iss <- lme4::lmer(legit_index_scaled_W ~ 
                                    black_W1 + age_W1 + education_W1 + female_W1 + income_W1 + 
                                    sc_info4_index_W1 + media_expo_index_W1 + 
                                    pol_trust_W1 + politicians_disagree_W1 +
                                    gm_family_gay_W1 + 
                                    WAVE*ideol_distance_W +
                                    (1|individualid), 
                                  data=aagmdatawaves,
                                  REML=T
)
summary(lmer6.idlxw.did.iss) 


# eof 